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EXECUTIVE SUMMARY 


The onset of laminar axisymmetric Rayleigh-Benard convection 
is investigated for a low-Prandtl number liquid metal in a 
cylindrical container. All surfaces are considered to be solid and 
no-slip. Two separate cases are examined for the thermal boundary 
conditions at the side wall, one with conducting and the other with 
insulated surface. 

The governing Boussinesq system is first perturbed and then 
simplified by introducing a Stokes stream function. Subsequently, 
a Chebyshev Galerkin spectral model is employed to reduce the 
simplified system to a system of first-order nonlinear ordinary 
differential equations. A local stability analysis determines the 
two values of the first critical Rayleigh number, Ra^^, for the 
insulated and conducting side walls. 

As expected, the conducting Ra^,]^ value of 2882.5 obtained from 
the present approach exceeded the corresponding insulated Ra^^ value 
of 2331.6. For the insulated case, an earlier study using a 
different numerical approach suggests that Ra^^ = 2261.9, while an 
experimental study measured Ra^^ = 2700. 
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CHAPTER 1. INTRODUCTION 


The homogeneity of semiconductor crystals is a primary 
concern of the electronics industry. The refinement of 
techniques for controlling the composition of semiconductor 
crystals during growth has been an integral part of recent 
advances in electronic microcircuitry. The difficulties and 
compromises inherent in growing homogeneous crystals have 
given impetus to both theoretical and experimental efforts to 
gain a better understanding of crystal growth phenomena. 
Recent theoretical work in the fields of nonlinear dynamics 
and stability theory has been promising in determining the 
onset of flow patterns in the liquid that affect the 
homogeneity of the resulting crystal. The computational 
difficulties of modeling a fluid on the route to turbulence 
remain formidable, however. 

Convection is the governing phenomenon in a number of 
crystal growth methods (e.g. , Ostrach, 1983). It is present 
in a wide variety of natural and industrial processes, among 
them some of interest to civil engineers. Thermal 
instabilities in the atmosphere (e.g., Lorenz, 1963) and in 
the borders of lakes (e.g., Horsch, 1988) are the product of 
convective forces. Secondary currents in stream meanders have 
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some similarities to the convection rolls present in the 
laminar regime of buoyancy-driven convection, the difference 
being one of exchange of momentum rather than heat. The study 
of convection in the more controlled conditions of crystal 
growth allows a better understanding of its effects in these 
other more complicated phenomena. 

A single semiconductor crystal can be formed from a 
liquid alloy, or melt, in the presence of a vertical 
differential temperature gradient. Although the most common 
configuration has a higher temperature at the top of the melt 
than at the bottom, the most interesting fluid phenomena are 
found in the "hot-on-bottom” configuration (e.g., Ostrach, 
1983) . The "hot-on-bottom" configuration is the case 
considered in the present study. 

The transfer of heat in the melt is controlled by 
conduction or by some form of convection, depending upon the 
strength of the temperature gradient. The physical properties 
of the crystal depend upon the type of heat transfer and the 
flow field in the melt during growth (Kim et al., 1972). The 
effect of convection on the crystal composition depends upon 
the degree of mixing in the melt (e.g., Muller et al., 1984). 
It is therefore possible to obtain some measure of control of 
the composition of the crystal by varying the flow pattern in 
the melt. Determining the particular flow regime in the melt 
at any time is a matter of determining the temperature 
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gradient at the onset of each of the different convective flow 
regimes . 

Single semiconductor crystals formed from alloys 
containing metals are a class of crystals with important 
industrial applications. Crystals such as lead-tin-telluride 
(LTT) have light-sensitive properties that have made possible 
the development of night -vis ion devices (Parker and Johnson, 
1981) , These crystals have unique problems associated with 
their growth. Because a melt containing a metal may be highly 
sensitive to changes in the temperature gradient 
(Krishnamurti, 1973) , control of the growth of the crystal can 
be problematic. In an effort to delineate the different flow 
regimes in such a melt, it is the purpose of the present study 
to determine the onset of the transition between the 
conductive and laminar convective regimes. 


1.1 The Growth of Crystals from the Melt 

There are a number of forces involved in the growth of a 
crystal . In a multi-component melt , convection may be induced 
by instabilities that are solutal as well as thermal (Crouch 
et al., 1985). In the absence of surface tension effects and 



4 


any horizontal temperature gradients, the flow in the melt is 
strictly some form of Rayleigh-BSnard (buoyancy-driven) 
convection (e.g.. Char Ison and Sani, 1970). 

Creating a multi-component crystal that is homogeneous 
can be an elusive undertaking. Although it is possible to 
grow a crystal in either a conductive or a convective flow 
regime, each heat transfer mechanism induces certain phenomena 
that create difficulties in controlling homogeneity. Under 
purely conductive heat transfer conditions, a multi-component 
crystal with a nearly constant concentration profile (Fig. 1) 
can be produced. The homogeneity of this crystal results in 
uniform electrical properties. However, large axial 
temperature gradients are necessary to avoid constitutional 
supercooling and dendritic growth at the crystal-melt 
interface (Tiller et al., 1953 and Mullins and Sekarka, 1964) . 
Sufficiently large axial temperature differences induce large 
radial temperature differences that are inherently 
destabilizing and initiate convection (e.g., Tritton, 1988). 
Convection that is unsteady produces a varying concentration 
profile in a multi-component melt (Fig. 1) . The more vigorous 
the convective mixing, the greater the occurrence of 
backmelting at the melt-crystal interface and the possibility 
of the formation of striations of differing compositions 
(e.g., Kim et al., 1972). Determination of the axial 
temperature gradients at which the flow regime transitions 
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occur is therefore essential for balancing the deleterious 
effects of these two heat transfer mechanisnis. 

In the present study the melt consists of the single 
component tin. In a single-component melt only thermal 
instabilities are significant, allowing this type of 
instability to be examined without the interference of most of 
these other complex phenomena. 

The specific response of the liquid melt to the forcing 
temperature gradient depends upon three parameters; (l) the 
strength and orientation of the temperature gradient, measured 
by the nondimensional Rayleigh number, Ra, (2) the relative 
rates of molecular diffusion of heat versus momentum by the 
fluid, measured by the Prandtl number of the fluid, Pr, and 
(3) the relative dimensions or aspect ratio 7 of the fluid 
container as well as its shape (e.g., Higgins, 1987). For a 
given fluid and container geometry, as the Rayleigh number is 
gradually increased, the melt passes from the conductive heat 
transfer state through a sequence of convective flow regimes; 
laminar, periodic, quasiper iodic, and turbulent (Krishnamurti, 
1973) . The Rayleigh numbers at which these transitions or 
Rayleigh instabilities occur are called critical Rayleigh 
numbers . 

The liquid metals used to grow semiconductor crystals are 
characterized by low Prandtl numbers, typically less than 0.5 
(e.g., Knuteson, 1989). Flow regime transitions for these 
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low-Prandtl number fluids occur over a narrower range of 
Rayleigh numbers compared to fluids of higher Prandtl number 
(Krishnamurti, 1973) (Fig. 2) . In the present study, a tin 
melt of Prandtl number 0.01 is heated from below in a closed 
cylindrical container with height equal to its radius. The 
first critical Rayleigh number marking the onset of the 
transition from conduction to the laminar convection regime 
for this set-up is determined. Although this critical 
Rayleigh number is independent of the Prandtl number and 
consequently of the fluid (e.g., Krishnamurti, 1973), the 
techniques used in this study illustrate the basic method in 
the determination of the entire series of critical Rayleigh 
number values for a buoyancy-driven flow. 


1.2 The Vertical Bridcfman Technique 

A common method for growing single crystals is the 
vertical Bridgman technique (e.g. , Carlson et al., 1984). In 
the ”top-down" vertical Bridgman configuration, a cylindrical 
ampoule of the liquid metal is moved upward through a linear 
vertical temperature gradient that is hotter on the bottom 
(Fig. 3) . The sides surrounding the ampoule are insulated, 
and the top and bottom surfaces of the ampoule are conducting. 
The crystal grows from the top of the ampoule down. The shape 
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of the crystal-melt interface is curved, with the fastest 
growth generally occurring at the sides. The Bridgman 
apparatus effectively makes the melt a closed system, so that 
there are no free surfaces and no destabilizing surface 
tension gradients . 

The enclosed nature of the Bridgman furnace setup (Fig. 
3) as well as the opacity of liquid metals do not permit the 
usual flow visualization techniques in the experimental 
apparatus. Instead, other techniques have been developed to 
assist in the determination of the temperature and velocity 
patterns for a liquid metal during a crystal growth 
experiment. For example, thermocouples arranged 
systematically about the cylindrical ampoule (Fig. 4) allow a 
temperature time series to be recorded (Fig. 5) . Deviations 
from the stationary conductive linear temperature gradient 
permit the large-scale temperature Variations to be detected. 
Because of the correlation between the temperature and 
velocity fields, especially in the vertical direction, an 
approximation to the large-scale flow field can be determined 
experimentally. This approximation is useful because it is 
precise enough to indicate the type of convective flow regime 
present. The Rayleigh instabilities may be found at the 
temperature gradients where there is a distinct and stable 
qualitative change in the temperature and flow fields. 
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1.3 The Solution Approach 

In the present study, the onset of laminar axisymmetric 
convective flow of liquid tin in a cylindrical container whose 
radius is equal to its height is investigated. The crystal- 
melt interface is assumed to be planar. The height of the 
melt coltimn is considered to be constant for the time scale of 
the study. The flow is considered to be axisymmetric, 
producing a rotational flow field in the shape of a torus. 
Because of the rigid surfaces of the container, the velocity 
is required to satisfy the no-slip condition — attain a zero 
value--at the solid boundary. For the thermal boundary 
conditions, the top and bottom are conducting faces. The side 
wall is insulated in the vertical Bridgman technique, but both 
insulated and conducting cases are considered in the present 
study as a check upon the validity of the model. Because of 
the loss of heat through a conducting side wall, the 
temperature gradient and critical Rayleigh number necessary to 
initiate convection should be higher than for the insulated 
case . 

In the present study, a Chebyshev-Galerkin spectral 
method is used to transform the governing nonlinear system of 
partial differential equations into a low-order nonlinear 
system of ordinary differential equations. A local stability 
analysis of this system is used to determine the first 
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critical Rayleigh number. 

In the present work, the problem is viewed first from a 
physical perspective. In Chapter 3, the notion of physical 
stability in Rayleigh-B^nard convection is explored, and the 
assumptions of the governing Boussinesq system are 
investigated. In Chapter 4, the Boussinesq system is 
perturbed and simplified. The dependent variables are 
represented by finite Chebyshev series, and a Galerkin 
spectral method is used to transform the Boussinesq system 
into an approximate system of first-order ordinary 
differential equations. 

In Chapter 5, the two types of instabilities of a 
mathematical system are examined. The nonlinear system is 
perturbed, linearized, and put into variational form. The 
conditions under which linear analysis of a nonlinear system 
is valid is presented. Such a linear analysis based upon the 
eigenvalues determines the critical Rayleigh number at the 
transition from conduction to laminar convective flow. Both 
insulated and conducting side walls are considered. 

The results are compared to a numerical study by Charlson 
and Sani (1970), who found the first critical Rayleigh number 
for the onset of laminar convection for both thermal boundary 
conditions for a range of aspect ratios. The results are also 
compared to the experimental results of Miiller et al. (1984) 
for liquid gallium. 
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Ficmre 1. Conductive (Diffusion-controlled) and convective 
(Completely-mixed) concentration profiles for a two-component 
melt (Knuteson, 1989) . 
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•fficpire 5. Temperature vs. time graph for periodic 
in a crystal melt. The TC-# indicates the 


flow field 
particular 


thermocouple (Knuteson, 1989) 



CHAPTER 2. LITERATURE REVIEW 


The difficulties that have been encountered in the study 
of Rayleigh-B^nard instabilities parallel those of two of the 
most intractable problems in fluid mechanics: hydrodynamic 
stability and the numerical modeling of turbulent flow. 
Despite the presence of only a single driving force, the 
nonlinear nature of Rayleigh“B4nard convection gives it a 
complex and often "chaotic" behavior that is difficult to 
describe and predict (e.g. , Berg§ et al., 1984). While the 
determination of instabilities in buoyancy-driven flows is 
intrinsically of interest because of the presence of Rayleigh- 
Benard convection in a number of important industrial and 
natural processes, the phenomenon has lately attracted more 
attention because the relative simplicity of the system lends 
itself to the study of the mechanisms involved in the 
transition to t\irbulence (e.g., Gleick, 1988). 

The development of mathematical tools for the study of 
Rayleigh-B4nard convection has not kept pace with the 
refinement of experimental efforts. Present numerical methods 
for solving the governing system of partial differential 
equations cannot adequately define the whole range of flow 
regimes and transitions (e.g., Tritton, 1988). In the field 
of nonlinear dynamics, methods for the study of mathematical 
stability have seen rapid growth in the past three decades but 
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are still undergoing development and refinement. 
Consequently, studies in Rayleigh-B4nard eonvection have been 
marked by the continued use of experimental studies to verify 
the theoretical approaches (e.g. , Tritton, 1988). 

B4nard (1900, 1901) was the first to conduct rigorous 
experiments on the problem of buoyancy-driven instabilities in 
a thin layer of spermaceti heated from below in a 
gravitational field. At the free surface of the fluid he 
observed a number of hexagonal convection cells with flow 
rising in the center of the cell and falling at the perimeter. 
It was found later by Block (1956) and Pearson (1958) that the 
cause of these cells was not vertical differences in buoyancy 
in the fluid layer, but inherently destabilizing surface 
tension gradients caused by surface temperature variations. 

From the theoretical standpoint, the first useful 
simplification of the equations governing buoyancy-induced 
flows in a thin fluid layer, the Boussinesq system, was not 
completely formulated until the turn of the century by 
Boussinesq (1903) . Boussinesq reasoned that for a small 
vertical temperature gradient, density variations were 
significant only in the buoyancy term of the Navier-Stokes 
equations. 

Lord Rayleigh (1916) assumed that the flow in B^nard's 
experiment was buoyancy-induced. He used a fluid layer thin 
enough to ignore the effect of the sidewalls on the flow 
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pattern. He performed a linear stability analysis on the 
Boussinesq system. Linear stability theory allows 
determination of the first critical Rayleigh number as well as 
the wavelength of the destabilizing perturbation, but it does 
not allow determination of the flow patterns of the velocity 
field (e.g., Knuteson, 1989). Subsequent investigations have 
refined considerably the linear stability of fluid layers of 
infinite lateral extent (e.g./ Chandrasekhar, 1961). 

Experimental work by Koschmieder (1966, 1967, 1969) 
showed that the side wall of the fluid container does indeed 
exert an inf luence on the convective pattern of the f low. For 
small vertical temperature gradients and certain values of the 
aspect ratio, the velocity field takes the form of convective 
roll cells in the shape of the surrounding side walls. In a 
cylindrical container under such conditions, he found 
axisymmetric rolls. Theoretical studies by Davis (1968) and 
Segel (1969) using large but confined geometries obtained the 
flow fields found in these experiments for rectangular box- 
shaped containers. 

The effect of the Prandtl number on flow regime 
transitions was investigated in exhaustive experimental 
studies by Krishnamurti (1970a, 1970b, 1973) . She found that 
flow regime transitions become increasingly more sensitive to 
variations in the Rayleigh number for fluids with Prandtl 
numbers approaching zero. In mercury, for example, with a 
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Prandtl number of 0.025, only turbulent flows are observed 
after the transition from conduction (e.g. , Tritton, 1988). 

Until the refinement of certain spectral methods of flow 
representation (Orszag, 1971a, b) , the theoretical study of 
Rayleigh-B€nard convection in cylindrical geometries was 
confined to the case of a side wall under "slip" conditions 
(e.g., Pellew and Southwell, 1940; Zierep, 1963; Ostrach and 
Pneuli, 1963; Liang et al., 1969). Zierep found that a no- 
slip side wall did not allow separation of variables. Orszag 
was able to develop flow representations using expressions 
other than Fourier series that allowed a rigid no-slip side 
wall and more realistic flow representation (Orszag, 1971a, b) . 
Charlson and Sani (1970) used a linearized Boussinesq system 
and a Rayleigh-Ritz method with Bessel function flow 
representation to find the first three critical Rayleigh 
numbers for axisymmetric flow. They considered both 
conducting and insulated side walls. Experimental work by 
Muller et al. (1984) with liquid gallium compared fairly well 
with their first critical Rayleigh number results. 

A straightforward numerical solution to the Boussinesq 
system of equations is possible using a number of techniques. 
There are two tasks that an efficient method must accomplish. 
The first is to efficiently model the whole range of flow 
regimes with high resolution on the route to turbulence. This 
flow modeling is complicated by the difficulty of representing 
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the dissimilar flow patterns and ever increasing modes of the 
flow on the route to turbulence. Tarman (1989) used two 
finite Fourier series of 106 and 297 terms to model turbulent 
Rayleigh-B4nard convection in a three-dimensional box of 
infinite horizontal extent and still was not able to 
accurately model vertical vorticity. Computational efficiency 
and rapid convergence are therefore paramount. Secondly, it 
is desirable to be able to detect flow regime transitions 
easily and precisely. Obtaining an efficient representation 
of the dependent variables temperature and velocity that can 
be varied easily with changes in other parameters such as the 
Rayleigh number to determine stability is difficult. 

Finite difference (e.g*, Muller et al., 1984), finite 
element (e.g., Carlson et al., 1985), and spectral methods 
(e.g., Canuto et al., 1988) are possible numerical approaches 
to Rayleigh-B§nard convection that can satisfy these 
conditions. Their application depends upon the particular 
characteristics of the problem under consideration. The 
presence and type of nonlinearities is often a deciding factor 
(e.g., Canuto et al., 1988). The primary difference between 
the finite difference and finite element methods versus the 
spectral methods is that the former divide the domain into a 
grid and solve the system for local grid points or elements, 
whereas the latter use global expressions to represent the 
dependent variables (e.g., Canuto et al., 1988). 
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To determine the instabilities, one approach is to mimic 
the sensitivity approach of the experimental methods. As the 
Rayleigh number is gradually increased, instabilities appear 
as distinct qualitative changes in the flow pattern solutions 
(e.g., Deane and Sirovich, 1990; Le QuerS and Alziary de 
Rochefort, 1986) . For computationally intensive models, this 
technique is tedious, as it necessitates implementation of the 
computational procedure with each new value of the Rayleigh 
number. This approach can also be complicated by the 
necessity to use long evolutions of the flow to avoid 
mistaking transient or intermittent flows for the true 
transitional flow. However, it may be the best approach to 
determining stability if the transformed systems of the 
spectral methods are difficult to analyze for stability. 

The use of Chebyshev functions to model flows began with 
Orszag (1971a, b) . He found it was possible to use Chebyshev 
series to model rigid side walls. In addition, he found that 
Chebyshev representations of dependent variables were suited 
to many hydrodynamic stability problems because of the 
accuracy and conciseness of Chebyshev approximations and their 
infinite-order accuracy compared to finite difference 
approximations. He used Fast Fourier Transforms to transfer 
between the physical and spectral spaces (1971b) . There have 
been a number of studies of eonvection to use Chebyshev series 
to model two-dimensional and three-dimensional flows. 
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Haldenwang (1986) found that Chebyshev spectral methods were 
suitable to modeling dissipation in a convective boundary 
layer at high Rayleigh numbers. Le QuerS and Alziary de 
Rochefort (1986) found the rapid convergence properties of 
Chebyshev functions to be advantageous in the problem of 
Rayleigh-Bdnard convection in a square cavity. Some 
difficulties that are encountered in these Chebyshev spectral 
flow representation techniques are the necessity of using an 
implicit or semi-implicit evaluation of the diffusive terms in 
the Navier-Stokes equation and the lack of natural boundary 
conditions for the pressure (Le Quer4 and Alziary de 
Rochefort, 1986). 

If the nonlinearities of the system are no more than 
quadratic, as they are in the Boussinesq system, then one 
efficient technique may be to use a Galerkin spectral method 
(Canuto et al., 1988). This approach transforms the 
Boussinesq system into a low-order system of ordinary 
differential equations. A stability analysis can then be used 
to determine the instabilities of this transformed system. 
The stability analysis of the dynamical system that results is 
a theory unto itself (e.g., Wiggins, 1990). The problem 
arises that a very complex nonlinear analysis may be 
necessary, one more difficult than the sensitivity approach of 
the other methods. 

For the present study, with its aim of determining only 
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the first transition to a relatively simple and approximately 
knovm laminar flow, a Galerkin method is well-suited to the 
problem. It is possible to represent the flow and the 
nonperiodic boundary conditions concisely, A local stability 
analysis is sufficient to determine the first critical 
Rayleigh number. The transformed system consists of only 
three equations, so that the eigenvalue equation is a cubic 
and the solution will be exact. 

The present study parallels the classic nonlinear dynamic 
analysis of Lorenz (1963) . Lorenz used a shorter form of the 
Fourier Galerkin spectral method representation of Saltzmann 
(1962) to analyze the stability of Rayleigh-Bdnard convection 
cells in the atmosphere. He considered periodic boundary 
conditions in a two-dimensional Cartesian framework and an 
infinite lateral extent. He discovered a rich behavior in the 
resulting low-order system of ordinary differential equations. 
His model was not extensive enough to accurately represent the 
whole range of actual flow regimes and Rayleigh instabilities 
(Nese, 1987) , but it gave impetus to the study of ''chaos,” the 
stability of nonlinear dynamical systems (e.g. , Gleick, 1987) . 

In the present work, the use of Chebyshev series in 
cylindrical coordinates to model the flow variables allows the 
imposition of a rigid, no-slip side wall in the cylindrical 
container. A Galerkin method is then used to transform the 
governing system of nonlinear partial differential equations 
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into a system of first-order nonlinear ordinary differential 
equations. A local stability analysis of the transformed 
system similar to that of Gelaro (1987) is used to determine 
the first critical Rayleigh number. 



CHAPTER 3. RAYLEISH-BENARD CONVECTION 


The growth of a crystal from the melt is a solidification 
process governed by the two heat transfer mechanisms in 
fluids, conduction and convection. The particular form of 
heat transfer depends upon the strength of the temperature 
gradient. As the temperature gradient is gradually raised 
from an initial zero value, the melt transfers heat by 
conduction and then by successively more vigorous forms of 
Rayleigh-Bdnard convection— laminar, periodic, quasiper iodic, 
and turbulent. A useful mathematical model must accurately 
represent both the physical complexity of each flow field as 
well as the onset of the transitions between them. 


3.1 The Physics of Ravleiah^BSnard Convection 

Rayleigh-B^nard convection is flow that is driven by 
buoyancy gradients in the fluid. As an illustration of the 
basic physical mechanism of this type of convection, consider 
a column of fluid heated from below in a gravitational field 
(Fig. 6) . The higher temperatures near the bottom cause the 
fluid there to expand, decreasing its density relative to the 
colder, more dense fluid above it. In this unstable 
equilibrium, the bottom fluid tends to rise and the top fluid 
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to fall relative to each other. Motion is initiated only if 
the driving temperature gradient is high enough to overcome 
the forces opposing motion. These forces are the frictional 
effect of the fluid viscosity and the damping effect of 
thermal conductivity, the tendency of the fluid to transfer 
heat by conduction rather than by convection. The direction 
of the flow rotation is determined by the initial conditions 
(Liang et al., 1969). If flow is initiated, the flow is 
horizontal at the top and bottom in order to satisfy 
continuity. The laminar convective flow pattern that is 
formed is a roll cell or BSnard cell. 

The simplest physical configuration illustrating the 
mechanisms of Rayleigh-Benard convection in the crystal melt 
consists of a layer of fluid between two extensive horizontal 
plates in a gravitational field (Fig. 7) . Consider that the 
temperature of the bottom plate T 2 is higher than the 
temperature of the top plate T^, T 2 > (Fig. 7) . For the 
conductive case the resulting temperature profile is linear. 
The ratio h/L « 1 so that there is no significant friction at 
the vertical boundaries and no influence of the sidewall on 
the flow pattern. 

Consider a temperature gradient that is marginally 
supercritical, sufficient to produce laminar convective flow 
from the conductive regime. Both conductive and convective 
heat transfer mechanisms are present in this flow, although 
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the latter is dominant. At the bottom of the fluid layer, 
heat from the bottom plate is transferred to the adjacent 
fluid by conduction. The increased buoyancy of this fluid, 
caused by the local decrease in its density, induces it to 
rise. Convective flow then transfers the heat to the vicinity 
of the top plate, where it is transferred from the fluid to 
the top plate again by conduction* The fluid, having lost 
much of its heat, falls and is heated on its way down. When 
it reaches the vicinity of the lower plate the whole process 
is repeated. 

The linear temperature gradient typical of conduction is 
modified by convection (Fig. 7) . The temperature in the 
convective region tends to a uniform temperature distribution, 
cooler fluid at the top being heated, warmer fluid at the 
bottom being cooled. The more turbulent the flow, the more 
uniform the temperature distribution tends to be in this 
region (e.g., Tarman, 1989). Temperature gradients in the 
conductive regions at the top and bottom boundaries are steep. 

For the marginally supercritical temperature gradients 
typical of laminar convective flow, the number of parallel 
rolls that are formed in such a thin fluid layer depends upon 
the aspect ratio of the fluid container. The vertical 
dimension of these single rolls is determined by the scale of 
the height of the fluid layer, h (Shirer, 1987) (Fig. 7) . 
This vertical dimensioning is valid even for containers with 
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h/L > 1 (e.g. , Tritton, 1988) (Fig. 6) . For a marginally 
supercritical temperature gradient, a single convection cell 
will be present in a cylinder with h/L > 1 (e.g. , Miiller et 
al. , 1984) . 

Although the flow field may become more complex with 
increasing temperature gradient, these basic heat transfer 
mechanisms of conduction at the boundaries and some form of 
convection in the middle remain valid. The flow fields for 
the other flow regimes are not so simple, but they can be 
thought of as superpositions of many different modes of this 
basic flow. 

3.2 The Boussinesa System 

The system of equations governing Rayleigh-B4nard 
convection consists of the continuity (conservation of mass) , 
conservation of energy, and Navier-Stokes (conservation of 
linear momentum) equations. In cylindrical coordinates, these 
equations are: 

|| + V-(pv) = 0 (3.1) 

= -(Vg) - - (t;Vv) 


(3.2) 
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(3.4f) 


where 


- heat capacity at constant volume, per unit mass 
- acceleration of gravity 
p = pressure 

q - energy flux relative to mass average velocity 
T = absolute temperature 
t = time 

V = velocity vector = (v^,Vq, v^) 


H - dynamic viscosity 

p = fluid density 

X - viscous stress tensor 

Tjj- = viscous stress tensor component 

( ) ^ = per unit volume 
( ! ) = scalar product of tensors 

This system is solved for the dependent variables v, p, and T 
in order to determine the velocity, pressure, and temperature 
fields. 

This system of three partial differential equations is 
nonlinear and coupled and is consequently difficult to solve 
numerically. The system can be simplified for buoyancy-driven 
convection under certain conditions to the Boussinesq system 
(Boussinesq, 1903) . The derivation is based on the assumption 
that for a small temperature gradient in the fluid layer, the 
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variation in the density is small but sufficient to cause flow 
due to buoyancy differences (e.g., Drazin and Reid, 1981). 
The density is assumed to be constant except in the buoyancy 
term pg^ of Navier-Stokes equation 3.3c, where the density is 
assumed to vary as a linear function of the temperature, 

Ap = -apgAT (3.5) 

where 

p = fluid density 
Pg = reference density 

a = coefficient of thermal expansion of the fluid 
At = temperature difference in the fluid layer 

other assumptions that are made in the derivation are 
(e.g., Knuteson , 198 9 ) : 

1. Gravity is the only external force and is 
directed downwards. 

2 . Accelerations in the fluid are small compared 
with g. 

3. All physical properties of the fluid are 
constant 

4. There is no internal heat generation. 

All these assumptions are satisfied in the system considered 
here. 

Rayleigh-B^nard convective flow in a low-Prandtl number 
crystal melt in the Bridgman configuration is therefore 
governed by the Boussinesq system with the following specific 
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conditions? a cylindrical geometry, no-slip conditions at all 
surfaces, conducting top and bottom, an insulated side wall, 
and a planar surface at the crystal interface. In three- 
dimensional cylindrical coordinates, the nondimens ional form 
of the Boussinesq equations (continuity, conservation of 
energy, and Navier— Stokes) is (e.g., Knuteson, 1989); 


V-v* = 0 

(3.6) 


(3.7) 

at* 


JL ( Jz! + v*-Vv*) = -Vp* + V^v* + Rar *§2 

(3.8) 


Pr at* 
where 

§ 2 , unit vector in the z-direction 

g, acceleration of gravity 

h, height of the cylinder 

p (r*, 0,z*,t*) , pressure 

Pr, Prandtl number = — 

K 

R, radius of cylinder 

r*, nondimensional radial coordinate 

Ra, Rayleigh number = - 

T* ( r * , 0 , z * , t * ) , temperature 

T 2 , temperature at the bottom boundary 

temperature at the top boundary 
AT, temperature difference in the fluid layer = T 2 - 
V* {X* ,6 , z* , t* ) , velocity vector 

z*, nondimensional vertical coordinate 

а, volumetric thermal expansion coefficient 

y, aspect ratio = h/R 

б , angular coordinate 
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K, thermal diffusivity 

p , kinematic viscosity 

~y/2^z*^yl2 cylindrical coordinate system 

For the nondimensionalization, velocities are scaled with 

K/h, temperature with AT, pressure with pK^fh^, time with 


h^/ic, and lengths with the height of 

the fluid column, h: 

JC 

(3.9a) 

r - r. 

(3.9b) 

ITT* — 2 

AT 

p* = 

PK^ 

(3.9c) 

II 

H 

(3.9d) 

_ txr 
h2 

(3.9e) 


(3.9f) 


For this particular problem, with the radius of the 
cylinder equal to its height, the aspect ratio y is one. The 
velocities are zero at the rigid, no-slip boundaries: 
v*(l,0,z*,t*) = v*(r*,9,-y/2,t*) = v*(r%^, 7 / 2 ,t*) 

= 0 (3.10) 

The top and bottom faces are conducting surfaces and are at 
constant temperature: 
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T* {X* ,9 C^elL, Ci>0»K (3.11a) 

T*(x*,e,-y/2,t*) = C 2 , C 2 €R, C2>0»K (3.11b) 

For the side wall, two cases are considered, insulated 

(l,0,z*,t*) =0 (3.12) 

dr* 

or conducting surfaces; 

T*(l,0,z*,t*)-C, C€R, C>0®ir (3.13) 


It is evident from equations 3. 6-3. 8 that dynamic 
similarity is attained by using the nondimensional parameters 
Ra and Pr. It is not possible to use the Reynolds number for 
dynamic similarity because of the absence of a characteristic 
velocity in the conductive regime. 

The Rayleigh number is a relative measure of the forcing 
mechanisms that drive Rayleigh-B^nard convection and the 
opposing mechanisms that damp the flow. The buoyancy 
differences in the fluid layer that create the thermal 
instability are the result of density variations. These 
variations are caused by the temperature difference AT and the 
consequent thermal expansion of the fluid, indicated by a. 
The height of the fluid layer h is included in the formulation 
because it is essentially the temperature gradient that is the 
forcing mechanism. Two effects inhibit convective flow, the 
viscosity p, the resistance of the fluid to shear deformation, 
and the thermal diffusivity k, the tendency of the thermal 
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conductivity of the fluid to transfer heat by conduction 
rather and thus postpone the onset of convection. 

The Prandtl number, a property of the fluid, indicates 
the relative molecular transport rates of momentum versus 
heat. It is evident from the Navier-Stokes equations (eqs. 
3. 3-3. 4) that a low Prandtl number significantly amplifies the 
role of the inertial nonlinearity (e.g., Knuteson, 1989). 
This effect is felt from the initiation of laminar convection. 

The aspect ratio 7 indicates the relative size of the 
container. It is one factor that determines the number and 
size of the convection cells in the container. For laminar 
convection in the cylinder of aspect ratio one, only one 
concentric roll is anticipated (Fig. 8 ) . 

The Boussinesq system has no direct analytical solution. 
Numerical solutions are approximations of the true solution, 
and for this reason they are not able to accurately describe 
or predict the full range of flow patterns of the fluid (e.g., 
Tritton, 1988) . In addition, the theory of stability of 
systems of partial differential equations is incomplete. 
Determination of the stability of the Boussinesq system must 
be approached indirectly, either by transforming to an 
approximate system of ordinary differential equations that is 
amenable to stability analysis or by a parametric study to 
detect qualitative changes in the flow or temperature fields. 

In the present study the former approach is taken. A 
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relatively simple flow is approximated using finite Chebyshev 
series. A Galerkin spectral method is used to transform the 
Boussinesq system into a first-order system of three ordinary 
differential equations. A local stability analysis is then 
performed to determine the first critical Rayleigh number. 
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Ficmre 6. Rayleigh-B^nard convection in a fluid column. 


T2 > Ti. 
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conductive temperature profile 
convective temperature profile 


Figure 7. Parallel laminar convection rolls in a closed 
rectangular box with slip side wails. The convective 
redistribution of heat is evident from the temperature 
gradients in the fluid layer (from Shirer, 1987) . 
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Ficnire 8. Axisymmetric convection roll in cylinder of aspect 
ratio one (from Miiller et al., 1984) . 


CHAPTER 4. THE CHEBYSHEV-GALERKIN SPECTRAL METHOD 

Galerlcin spectral methods transform a system of partial 
differential equations into a low-order system of ordinary 
differential equations, a form more amenable to solution as 
well as to stability analysis. A Chebyshev-Galerkin spectral 
method is especially useful for modeling nonperiodic laminar 
flows in simple geometries (Orszag, 1971b) . Many such flows 
can be approximated by continuous polynomial functions, and 
because Chebyshev functions are essentially polynomials, 
dependent variables can be represented by finite series of 
Chebyshev functions. A Galerkin method employing Fourier 
series, although much easier to apply in such cases, would not 
allow accurate modeling of the nonperiodic flows and boundary 
conditions. The set of Chebyshev functions is also a member 
of a class of functions having the orthogonality properties 
necessary to implement the Galerkin method. 

The laminar axisymmetric flow of a low-Prandtl number 
liquid metal in a rigid, closed cylindrical container of equal 
height and radius is a case of a relatively simple nonperiodic 
flow with nonperiodic boundary conditions. The convection 
cell has a toroidal shape that is assumed to be symmetric 
about the horizontal mid-plane (Fig. 9) . Although the 
convergence of flows at the center of the cylinder displaces 
the rotational center of the cell outward to a point past R/2 
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(Muller et al., 1984)), both velocity and temperature with 
their boundary conditions can be represented fairly accurately 
by continuous polynomials and therefore by Chebyshev series. 

The physical transition of the fluid from the conductive 
to the convective regime cam be approached mathematically as 
a perturbation problem. The laminar convection state is 
treated as a perturbation of the basic conductive state. 
Analysis of the growth of the convective perturbation 
indicates the critical point of transition between the two 
regimes . 

The axisymmetric Boussinesg system is first perturbed. 
It is then simplified by transforming the Navier-Stokes 
equation into a vorticity equation and introducing a Stokes 
stream function to reduce the ntamber of dependent variables. 
A Chebyshev-Galerkin method is then applied to the resulting 
system of equations. Dependent variable Chebyshev series 
consisting of terms with separated variable expressions for 
time and the two spatial directions are substituted into the 
Boussinesg system. An ordinary differential equation of first 
order for each of the time-varying coefficients is obtained by 
multiplying a partial differential equation by the orthogonal 
test function of the desired time-dependent coefficient and 
integrating over the domain. Once the partial differential 
Boussinesg system has been transformed into a low-order system 
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of ordinary differential equations, the stability analysis can 
be performed. 


4.1 Perturbation of the Boussinesa System 


Perturbation of the Boussinesq system is useful because 
it allows the system to be expressed in terms of the 
convective flow variables. Normally, a perturbation analysis 
for a strongly nonlinear equation like the Navier-stokes 
equation is not straightforward. However, in this case, the 
zero velocity field of the conductive flow field allows 
nonlinear terms to be retained in the first order equations. 

The nondimensional cylindrical Boussinesq system from 
Chapter 3, 

V'V* = 0 (4-1) 


=s <p* 

dt* 


(4.2) 


_1_ + v**Vir*) = -Vp* + V^v* + RaT*% 

Pr dt* 


(4.3) 


is used here. Assuming axisymmetric flow, the above system of 
equations can be simplified by considering the following 
properties of this flow, which hold by definition: 
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v*^(l*,0rZ*,t*) - 0 (4.4) 

The first step in the perturbation procedure is to 
perturb the dependent variables velocity, pressure, and 
temperature. These variables are expressed as the sum of a 
basic state, conduction (O-subscript terms), and a small 
convective perturbation of this basic state (primed terms) : 


v*(r*,z*,t*) 

= v$ + v*'(r*,z*,t*) 

(4.6) 

p*(r*,z*,t*) 

= Po(z*) + p*r(r*,z*,t*) 

(4.7) 

r*(r*,z*,t*) 

= r$(z*) + r*'(r*,z*,t*) 

(4.8) 


In the conductive state, the velocity vector is zero, the 
temperature variation is linear, and the pressure gradient is 
hydrostatic: 

vj = 0 (4.9) 

Pa(z*) = P2 “ fg h3/«2) (z* + 1/2) (4.10) 

To(Z*) = T2 - (Z* + 1/2) (4.11) 

The terms P 2 and T^ are the pressure and temperature, 
respectively, at the bottom of the cylinder. 

Convective flow occurs if the perturbations are not 
damped over time. In that case, the perturbation variables 
represent the convective flow field. 

The expressions for v*, p* , and T* are substituted into 
the Boussinesq system. After some algebraic manipulations and 
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deletion of primes and asterisks, the system for the perturbed 
variables becomes; 

V-v « 0 (4.12) 

+ vVr - Vg = (4.13) 

+ vVvj » -Vpo - Vp + V^v + (RaTo + RaT) §g (4.14) 

Note the retention of the nonlinearities in the Navier- 
Stokes and energy equations. In this perturbed Boussinesq 
system as in the unperturbed, it is again apparent that a low 
Prandtl number exerts a significant nonlinear effect due to 
the inertial term 

(4.15) 

Pr 

This effect becomes important as soon as convection is 
initiated. In any flow regime with nonzero finite velocity 
vectors, the Navier-Stokes equations for low-Prandtl number 
fluids are highly nonlinear (e.g., Knuteson, 1989). 


4.2 Simplification of the Perturbed Bousainesa System 

The stability analysis of the perturbed Boussinesq system 
involves the solution of a variational equation that is a 
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polynomial. The degree of this equation can be decreased and 
its solution facilitated by reducing the number of dependent 
variables in the perturbed Boussinesq system. In addition to 
decreasing the size of the resulting system of ordinairy 
differential equations, this simplification also eliminates 
the need to determine the pressure. 

Both pressure gradient terms -VpQ and -Vp' as well as the 
conductive buoyancy term RaToS^ in the perturbed Navier-Stokes 
equation can be eliminated by taking the curl of the equation 
and thereby transforming it into a vorticity equation: 


X - V X (V^v + Ra-TB^) (4.16) 

A nondimens ional Stokes stream function ^ is introduced 
to reduce the number of dependent variables from three 
(Vr/Vg/T) to two (^,T) . The ^ is made nondimensional by 
scaling it with the thermal diffusivity, k, and the height of 
the cylinder, h. The nondimensional radial and vertical 
velocities Vj. and v^ may then be expressed in terms of \pi 




1 difr 
r dz 


(4.17) 


^2 = 


1 dt 
r dx 


(4.18) 


Because this stream function satisfies the continuity equation 
exactly. 
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-Ra 

\ r 5z r dr 


8Ta 


in which the Jacobian J(fi,f 2 ) is 


JCfi/fa) = 


dr dz dz dr 


(4.21) 


(4.22) 


This system of two equations in ^ and T requires ten 
boundary conditions for the problem to be well posed. In the 
energy equation (eq, 4.20) , which is second-order in T in r 
and z, four boundary conditions are necessary to uniquely 
determine the coefficients of the particular solution for T, 
two from the r-variable and two from the z-variable. In the 
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same manner, six boundary conditions are needed for the 
solution of ^ in the third-order vorticity equation (eq. 4.21) 
in r and z. 

The ten boundary conditions for this problem are based on 
the model of the cylindrical ampoule used in the Bridgman 
method as well as one characteristic of axisymmetric flow. 
The ampoule is considered to have rigid, no-slip surfaces: 


-f -f .0 


(4.23a) 


^(l,z,t) =-|^(r,— |,t) =-|^(r,-|,t) = 0 (4.24b) 

oz oz 2 az 2 

The radial velocity on the center line of the ampoule is zero: 

^{0,z,t) = 0 (4.24) 

az 

The top and bottom are conducting surfaces. They conduct 
heat into and out of the fluid at the fluid-surface interface, 
but they maintain a constant temperature. Mathematically, the 
temperature perturbation is zero: 

T(r,-|,t) =r(r,l,t) =0 (4.25) 

For the thermal boundary conditions at the side wall, two 
cases are considered, in the first case a side wall that is 
conducting. 


r(i,z,t) = 0 


(4.26) 
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and in the second, a side wall that is insulated and allows no 
transfer of heat across the interface. 


dT 

8r 


(l,z,t) = 0 


(4.27) 


4.3 Chebvsfaev Series Representation of Dependent Variables 

In the Chebyshev^Galerkin spectral method, the dependent 
variables are represented by finite Chebyshev series, A 
separation of variables approach is used. Each term in the 
series is the product of two spatial ly-dependent expressions 
and a time-dependent amplitude coefficient. Successive terms 
represent modes of the flow, and the series is truncated when 
the desired number of modes is represented. The stream 
function and the temperature T in the perturbed Boussinesq 
system may be represented in the following manner; 

i/r{x, z,t) ^j (z) i, je{0, 1,2, . . . } (4.28a) 


TlTij ( t ) [aiTi (r ) ] [ bjT j ( z ) ] a^ , eR (4.28b) 
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T(x,Z,t) = ?? 1^2 (t) Tj^Cr) Tj (z) lC/ l6{0/X/ 2/ • ■ • } (4 •29a) 

= EE tCkTk(r)] [diTi(z)] Ci,,di€a (4.29b) 

k 1 

in which the Chebyshev function T„ is defined by: 

Tjj(x) = cos (n arccos x) x€[-l,l] , ne{0,l,2, . . .} (4.30) 

The stability of a single mode, laminar axisymmetric 
flow, is being considered in the present study. The relative 
simplicity of this flow allows the number of terms and the 
specific coefficients in these Chebyshev series to be 
determined fairly precisely. This axisymmetric flow can be 
approximated easily by using Chebyshev functions, to the 
extent that the spatial coefficients aj^, bj, c^, and dj^ can be 
determined. Even though a large number of terms is required 
to represent this flow, the analysis is simplified by the fact 
that the presence of a single mode requires that all 
constituent ” submodes” move together. The same temporal 
coefficient must then be common to all terms. The most 
significant result of this conciseness of the flow 
representation is that the stability analysis is greatly 
simplified. 

By way of contrast, a model that investigated the series 
of flow regime transitions on the route to turbulence would 
have to incorporate a very large number of terms with 
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differing temporal coefficients. The coefficients bj, Cj^, 
and dj^ would not be able to be specified and would be 
incorporated into the appropriate temporal coefficients. Only 
in this manner could the model represent the multiplicity of 
flow regimes as well as account for the increasing nximber of 
superimposed modes within each flow regime as it becomes more 
complex. 

The spatial distribution of ^ in the r and z planes is 
determined by the continuity requirement V-v = 0 as well as by 
the boundary conditions for Vj. and v^., obtained from the 
boundary conditions for With these constraints, it is 
possible to represent the entire flow field by a single term 
for \j/ in the stream function series. 

Letting ^ be the expansion 

= ^^(t) ^(i) fiz) , (4.31) 

then by using the Stokes stream function definitions, the 
distributions of Vj. and in r and z can be represented in 
terms of their corresponding ^ distributions: 

vv (r,z,t) = (4.32a) 
r oz 


= -^(t) 



difriz) • 
dz 


(4.32b) 


= Vj.it) Vj(r) Vjiz) 


(4.32c) 
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_ 1 d4r 
r 8r 


(4.33a) 

= ^(t) 


(4.33b) 

= v,(t) 

vjr ) v^ iz) 

(4.33c) 


It is possible to approximate the velocity distributions by 
using polynomials. Use of these direct relations between the 
velocities and ^ then allows determination of approximate \j/ 
spatial terms in Chebyshev series. 

For example, examination of the axisymmetric flow field 
(Fig. 10) reveals that the radial variation in v^{r ,z ,-t) for 
re [0,1] may be approximated by a cubic polynomial (Fig. 11). 
This rotational flow is assumed to move downward at the center 
of the cylinder. This direction of flow is not critical. 
Four conditions that specify this cubic curve include the 


velocity boundary conditions: 

• = 0 (no-slip condition at side wall) (4.34) 

• Vp(l,z,t) = 0 (no-penetration condition) (4.35) 

• 9v_/3r (0, z,t) = 0 (maximum v^ value at center 

line) (4.36) 

• v^{0,z,t) = -1 (arbitrary intercept) (4.37) 


The derivation of Vj.(r) depends directly upon v^(r) , both 
being functions of The Vj. condition above (eg. 4.35) is 
necessary to ensure that the form of v^ allows this condition 
to be met upon derivation of Vj.. 
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Letting v^(r) be the following general cubic polynomial, 

Vj,(r) = ar^ + br^ + cr + d, (4.38) 

the intercept condition (eq. 4.37) determines d = -1, and the 
maximum condition (eq. 4.36) determines c — 0. Substituting 
these values, the no-slip condition (eq. 4,34) becomes 

a+b = l (4.39) 

For the no-penetration boundary condition, v^(r) must first be 
derived, employing the mutual dependence on yf/. The variation 
of ^ in r is first found by taking the polynomial for V 2 .(r) 
and Integrating to solve for i^(r). 


v,(r) = 

1 a^(r) 
r 3r 


(4.40) 

iHr) -j 

frv^(x) dr 


(4.41a) 

= 

|*r(ar^ + br^ - 

l)dr 

(4.41b) 


^7- 5 + t 4 _ 

5 4 

-i-r^ + e 
2 

(4.41c) 


The variation in r for v^ then becomes (Fig. 12) : 
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Vj(r) = Y (4.42a) 

= ■*■ (4.42b) 

4 iu X 

The constant of integration "e" must be zero in order for 
Vj.(r) to be defined at r = 0. The no-penetration condition 
then becomes: 


a 

5 



(4.43) 


Solving for the coefficients a and b using the above equation 
4.43 and equation 4.39, a - -5 and b = 6. The equation for 
v^(r) becomes (Fig. 11): 

Vj,(r) = -5r^ + 6r^ - 1 (4.44) 

This curve models the salient physical features of the 
flow. The change in sign of the vertical velocity in this 
closed container indicates flow rotation. The location of the 
r-intercept to the right of r R/2 as well as the relative 
magnitudes of the maximum and minimum on the interval show the 
flow convergence at the center. The flow profile at the wall 
takes the approximate shape of a boundary layer. 

The equation for Vj.(r) is (Fig. 12) : 
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(4.45) 


The radial velocity v^ is zero at r = 0, and the profile near 
the center is similar to that of a boundary layer. 

The equation for ^(r) is: 


^(r) =* -r® + - —x^ 

2 2 


(4.46) 


For the purposes of the Galerkin method, V'(r) may be expressed 
in Chebyshev form as: 




^T2(2r-l) -^t,(2x- 1) -^To(2r-l)j (4.47) 

Expressing \p in terms of the quantity 2r-l for re [0,1] 
effectively translates the expression into the domain of all 
Chebyshev functions, [-1,1], and allows the orthogonality of 
Chebyshev functions to be used to simplify the integration of 
the equations in a later step. 

In a manner similar to the procedure above, the ^(z) , 
Vj.(z) (Fig. 13), and Vj,(z) (Fig. 14) terms for ze [-%,%] can be 
determined. Using the relationships in equations 4.32b,c and 
4.33b,c the equations for Vj.(z) and v^{z) can be related 
through \J/: 
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(4.48a) 


3y^(z) 

dz 


(4.48b) 


The following conditions are used to derive a quartic 
distribution for Vj^(z) and a cubic distribution for Vj.(z): 

• v^(r ,±h,t) — 0 (no-penetration) (4.49) 

• Vj.(r,±%/t) = 0 (no-slip at top and bottom boundaries) 

(4.50) 

• dv^/dz(r,o,t) - 0 (maximum value at midplane) (4.51) 

• v^(r,o,t) = 1 (arbitrary intercept) (4.52) 

• 3v3./3z(r,±.29,t) = 0 (maximum Vj. value) (4.53) 

The last condition for fixing the height of the maximum 
radial velocity Vj.(z) allows the thickness of the boundary 
layers at the top and bottom to be modified. The values of 
z=±.29 were chosen to smooth the transition of the maximum 
v^(z) value at 0.8R near the side wall (Fig. 11) to these 
maximum Vj.(z) values at z=±0.21h. The ±.29 values also allow 
integer values of the coefficients. 

After some algebra, the following distributions are 


derived: 
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fiz) - 16 - 8z^ + 1 


(4.54a) 


• 1 
128 


T4(2z) 


7 

32 


T2<2z) 


+ 


99 

128 


To(2z) 


(4.54b) 


v_ (z ) = 64 z ^ - 16 z 


(4.55) 


VjCz) * 16z^ - 8z^ + 1 (4.56) 

Note the antisymmetry of the v^{z) curve, the flow moving in 
towards the center at the top and outward at the bottom. 
Combining the two spatial expressions, \p becomes: 


[-- ^“T5(2r-1) 


— T 
128 * 


(2r-l) 




^To(2r-l)] [^T.(2 z) -^T,(2z) . 


99 

128 


To ( 2 z)] 


(4.57) 


in Which the subscripts i and j of denote the highest 
degree Chebyshev function in r and z, respectively. 

It may appear that the above Chebyshev series for t/' (eq. 
4.57) does not follow the form given in eqn 4.28b for the 
Galerkin method, in that there is only one temporal 
coefficient, ^ 54 (t) , for all eighteen distributed products 
T£(r)Tj(z). However, there is only one mode being considered 
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here, so that each constituent part of the flow moves together 
with the whole. The temporal coefficient is the same for each 
term and may be factored out. 

The flow field appears as in Fig. 15 for constant unit 
temporal coef f ic lent s . 

In the case of the temperature expansion, the choice of 
terms is based upon physical arguments from experimental 
findings. The two phenomena that are significant are the 
vertical redistribution of heat due to convection and the 
strong correlation between the temperature distribution and 
the vertical velocity (Shirer, 1987). The temperature 
function can be decomposed into a two-term expansion 
incorporating these features: 

r(r,z,t) = r^(r,z,t) +jr^(r,z,t) (4.58) 

The strong correlation between the temperature 
distribution and the vertical velocity allows use of a 
temperature profile similar to that of the vertical velocity 
v^. Two separate thermal boundary conditions are considered 
for the side wall. In the case of a conducting side wall, the 
temperature perturbation is zero at the wall, and the 
temperature distribution r^(r,z,t) (Figs. 16, 17) is assumed 
to be the same as that of v^ (Figs. 11, 14): 



57 


ri(r,z,t) = rMt) r^(r) r^Cz) (4.59a) 

= l34(t) [-5x3 + 6J.2 _ 2] [X6z* - 8z3 + 1] (4.59b) 

For the case of insulated sides, in which the derivative 
of the temperature with respect to r is equal to zero, a fifth 
degree polynomial that approximates the r-variation but with 
a zero r-derivative at r - 1 is chosen (Figs. 18, 19); 

r3(r,z,t) = T3(t) T^ix) T^iz) (4.60a) 

= r54(t) [6.0x5 - 11. 3x^ + 2.0x3 + 4.6x3 - 1] * 

[16 z^ - 8 z 3 + 1] (4.60b) 

The coefficients for this polynomial are rounded to the 
nearest tenth. In the attempt to construct a curve similar to 
the conducting case profile at every point except near the 
side wall, the following conditions were used; 


' r(0,z,t) = -1 (same T-intercept) (4.61) 

• T(.558,z,t) « 0 (same r-intercept) (4.62) 

’ T(l,z,t) = .28 (same maximum T value) (4.63) 

• r'(0,z,t) =» 0 (minimum r at r = 0) (4.64) 


• 3''(l,z,t) = 0 (r-derivative of T = 0 at r=l) (4.65) 

• r'(.528,z,t) = 2.02 (same slope at r-intercept) (4.66) 
The relative maximum over re [0,1] of the resulting curve is 
not at the side wall, but at r « .857. Its value of -.297 
does not vary sufficiently from the conducting side wall 
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maxlmiim of .280 to cause undue concern. The two different 
distributions are compared in Fig. 20. 

The second term of the temperature expansion represents 
the vertical redistribution of heat due to convection. 
Relative to the linear conductive temperature distribution, 
convection heats the fluid at the top of the container and 
cools the fluid at the bottom. This antisymmetric 
distribution is independent of r, and it can be represented by 
a cubic polynomial (Fig. 21) : 

r^(r,z,t) = ToaCt) [-32z^ + 8z] (4.67 a) 

= ro3(t) I-T3(2 z) + Ti(2z)] (4.67b) 

The final Chebyshev series for the two temperature 
expansions are, for conducting sides. 


r(r,z,t) - r 34 (t ) 


XT3(2r-l) - -^T2(2r-1) + -||Ti(2r-l) 


16 


T(2r-1) 


^T^i2z) - -|t2(2z) + ■|To(2z) 


To3(t) [-T3(2 z) +Ti(2z)] (4.68) 

and for the insulated side wall, 

r(r,z,t) = T^^it) [6. Or® - 11. 3r^ + 2.0r® + - 10]* 

[16z^ - 8z^ + 1] + To^it) [-32Z® + 8z] (4.69a) 
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= 1^4 ( t) 



(2r-l) 


+ 


37 


1280 


T4 


(2r-l) - 


3725 

32000 


T3 


(2r-l) 


_ 925 „ 
8000 ^ 


(2r-l) 


483 

640 


T, (2r-l ) - 


11857 rr, 
1280 ° 


(2r-l) 


■-|t.(22) 


T,(2z) + 4 t.(22) 


2 2 


8 


ro3(t) [-T3(2 z) +Ti(2z)] 


(4.69b) 


4.4 Chebvshev-Galerlcin Spectral Method 

The Galerkin spectral method converts the perturbed 
Boussinesq system into a system of ordinary differential 
equations in the temporal coefficients. The reason for this 
transformation is that stability theory for systems of 
ordinary differential equations is much easier and more 
refined than for systems of partial differential equations. 
Although this system models the flow in the conductive and 
laminar flow regimes, it is the mathematical stability of the 
system and the physical stability of the modeled flow itself 
that are actually being considered. 

The inaccuracy inherent in the representation of the 
dependent variables \p and T by approximate Chebyshev series is 
mitigated by the weighted residual technique of this method. 
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The general procedure (e.g., Canuto et al., 1990) is as 
follows. 

Let the perturbed Boussinesq system be represented by 
- M(u) = 0 (4.70) 

in which 

u(r,z,t) = (4.71) 

The Chebyshev series approximation for u is 

u"(r,z,t) = 3jt(t) Ujt(r) Ujfc(z) (4,72a) 

k-l 


N 

’ p 

' R 


EbiTid) 

£ciTi(z) 



b-o J 


(4.72b) 


These expansions may not satisfy eg. 4.70 exactly over the 
whole interval, the residuals 

-^-M(u") (4.73) 

not being everywhere zero. The optimal solution is obtained 
by requiring that the integrals of the weighted residuals. 



- M(u^ 


^(r)^(z) didz 


0 , 


k = 1, 


./N 


(4.74) 


vanish over the interval. 
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The residual is weighted by the test functions 0(r) and 
0 ( z ) . These two functions are chosen to be orthogonal to the 
corresponding spatial expansions U]^ over the interval; 


f^Ujt(r)0i(r)dr = Ci€]R,Ci>=0 

Jo 

.1 

J 2 ^u*(z) 0 i(z) dz = 026^1, CjeM/ Gj’^O 

"■2 


( 4 . 75 ) 


( 4 . 76 ) 


in which 5]^^^ is the Kronecker delta. 

In the Galerkin method, test functions are chosen to be 
the same type of functions as in the spatial expansions, in 
this case, Chebyshev functions. For a single Chebyshev 
function Tj^(x) , the test function is 


Tj (x) 


( 4 . 77 ) 


The weighted orthogonality condition is: 



Tj(x) 



0 

JL 

2 

n 


i’^j 
i=j#0 
i=j =0 


( 4 . 78 ) 


The weight is the radical denominator. For a spatial 
expansion consisting of a finite Chebyshev series, such as 


U(X) = 


i»0 


( 4 . 79 ) 


one possible test function is: 
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0(X) = 

^l-x^ 


(4.80) 


The result of the orthogonality in this case is 


P u (x) 0 (x) dx = 

•l 

52 biTi (x) 


J-i J 

-1 


■/1-X2 


(4.81a) 


= It 


2.1 /• l .2 

0+-— (bi 


hi 


+ 



(4.81b) 


This particular test function is chosen for two reasons. The 
first is that having the sum of the squares of the Chebyshev 
coefficients b^ (eg. 4.81b) ensures that the entire set of 
temporal coefficients aj^(t) of equation 4.72 will always 
appear in the equation, never taking a value of zero. The 
second reason is that each constituent part is certain to be 
represented. This type of test function is applied 
consistently to all three derived equations. 

Once these test functions have been determined, the 
system of integral equations can be solved, resulting in a 
low-order system of ordinary equations in the temporal 
coef f ic lent s : 



dt 


- f (a;,) = 0 
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k = l,2, ...,N (4.82) 

This Chebyshev-Galerkin spectral method is applied first 
to the case of conducting side walls. This situation is 
algebraically easier than the insulated side wall case because 
of the less complex form of the temperature expansion. The 
T 34 (t) equation is found first to illustrate the method. 

The T 34 (t) equation is derived from the perturbed energy 
equation, which contains the time derivative of T and 
therefore the time derivative of its component r34(t) . 

The expansions for the stream function ^ and the 
temperature for conducting side walls, 

(t) ?^(r) ?^r(z) (4.83a) 




-^T5(2r-1) - + 


^T3(2r-1) +^T3(2r-l) - .^T3(2r-1) 


128 


To(2r-l) 


^T,(2z) - ^T3<2 z) . 


99 

128 


Tq(2z) 


(4.83b) 


and 
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r(r, 2 /t) * r 34 (t)rHr)r^( 2 ) +ro3(t)r2(2) (4.84a) 

= f-^T3(2r-l) -^T2(2r-1) + ||T^<2r-l) - 

■^T(2r-1)‘ 'iT4(22) - -|t2(22) + 1To(22)J + 

ro3(t) [-T3(22) +Ti( 22)] (4.84b) 

are substituted into the perturbed energy equation, 

4? + — J(^, r) -iM -V2r= o (4.85) 

at r r or 

The equation is then multiplied by the test functions for 
2*34 (h) and integrated over the domain: 

nA'i? * ^]^<’'>*<^>}drdz = 0 (■‘•86) 

~~2 ° 

where the test functions for r 34 (t) are 



(4.87b) 


and 
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0(Z) = 


2ri<z) 

Vl-( 2 z )2 


(4.88a) 


-iT^Caz) - -iT2(2z) + |T<,(2 z) 
Vl-(2z) 2 


(4.88b) 


The coefficient ”2" in the numerator of the test functions is 
the constant resulting from the transformation of eq. 4.78 in 
X to one in either 2r-l or 2z. 

The software package "Mathematica'' was used to perform 
the integration. The capabilities of this software are such 
that each term must be evaluated separately. The resulting 
equation for the time variation in r 34 (t) is: 

TLit) = - ^409i^ ^4(t) - (4.89) 

Note that the equation is nonlinear, and that all terms have 
been retained from the original equation. 

This procedure is repeated to obtain the ro 3 (t) equation 
from the energy equation and the ^(t) equation from the 
Navier-Stokes equations to complete the system: 

ro3(t) =-ro3(t) * (4.90) 


^«(t) = + i|prRar34(t) (4.91) 

In the T'o 3 (t) equation, note that the integration of the 
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(-1/r) (d^/dr) term produces a zero. In the (t) equation, 
the term derived from the vVv term of the perturbed 
Boussinesq system is zero. The absence of this nonlinear term 
in the case of a low-Prandtl number fluid would be significant 
if the stability of any convective flow regime with its finite 
velocity vectors were being considered. In the present work, 
however, conductive stability with its zero velocity vector is 
of interest, and the loss of this inertial term has no effect. 

These same features can be seen in the case of insulated 
side walls. The same Chebyshev-Galerkin procedure is 
followed, the difference being the use of the insulated form 
of the temperature expansion and its test function. 


ilr54it) = Pr ^54 (t) + Pr Ra 


(4.92) 


TUit) = 


114636 

1205993 


^54 - 


38915152.,, 

1361605 


3209800 

6029965 


^54(t)ro3{fc) 


(4.93) 


r'aCt) = -48ro3(t) + mil ^5,(t) (4.94) 

In form, this system is similar to that of the conducting side 
wall case, the only difference being the replacement of the 
Tq 2 ('^) term in the nonlinear term in the equation 

(4.90) with T 54 (t). The most noticeable difference is that 
the magnitudes of the coefficients of T^^(t) are relatively 
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larger in this insulated set of eguations than those of the 
coefficients of r34(t) in the conducting case. This increase 
may introduce more damping of the growth of the stream 
function variable and consequently of the flow. 

Both of these systems of first order ordinary 
differential eguations are now in a form amenable to stability 
analysis. The small number of eguations is a reflection of 
the relative simplicity of the laminar axisymmetric flow. In 
the stability analysis of the next chapter, the critical point 
of transition from the conductive regime to the convective 
regime is determined. 
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Figure 9. Toroidal shape of axisymmetric convection cell. 





RADIAL VELOCITY, Vr(r) 


7 



T?-i miT*« 1 nn n*f vr . 


r 



VERTICAL DISTANCE FROM CENTER, 


73 



Picmre 14. z-variation of 
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Fxcrore IS. Axisyinmetric flow field for ^(r,z) in half -section 
of cylinder with 7-I. Flow rotates counterclockwise, falling 
at the center. 





VERTICAL DISTANCE FROM CENTER, 


76 



Figure 17. z-varlation of conducting T^(r,z,t) expression. 
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Ficmre 18. r-variation of insulated T^(r,z,t) expression 
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Ficmre 20. r-variations of insulated, conducting T^(r,z,t) 
expressions. 




CHAPTER 5. LOCAL STABILITY ANALYSIS 


The crystal melt exhibits a niamber of different flow 
patterns on the route to turbulence, from the conductive state 
of quiescent flow, through laminar and periodic convective 
flows, to full-scale turbulent convection (e.g., Knuteson, 
1989). The point of transition at which the onset of a new 
flow regime occurs is called a bifurcation. The Rayleigh 
mamber values at these points of bifurcation are known as 
critical Rayleigh numbers, and they are numbered in the order 
in which they occur. In this study, the critical Rayleigh 
number at the onset of laminar convection, the first critical 
Rayleigh number, Racj^, is determined. 

Each of the different flow regimes is "stable" over its 
particular range of the Rayleigh number. A stable flow regime 
is one in which the physical perturbations induced by changes 
in the forcing temperature gradient are sufficiently damped by 
viscosity and thermal conductivity to maintain the distinct 
qualitative characteristics of the flow. Flow instability at 
the bifurcation is caused by the undamped growth of 
perturbations, resulting in time in a qualitatively different 
flow pattern (e.g., Gelaro, 1987). To determine the stability 
of a flow or the stability of its mathematical model is to 
determine the bifurcations. 

The flow patterns in the melt are characterized by 
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greater complexity and loss of regularity as the Rayleigh 
nxunber is increased. The stability of the corresponding 
solutions of the mathematical model can themselves be 
categorized by their degree of complexity. Bifurcations are 
either local or global, depending upon the degree of 
regularity of the temporal behavior of the solution (e.g. , 
Thompson and Stewart . , 1986) . For the two systems of ordinary 
differential equations derived in the last chapter, the local 
bifurcation at the transition from conduction to laminar 
convection will be examined. The stability analysis is then 
said to be a local stability analysis. 

The local stability analysis of systems of ordinary 
differential equations that are nonlinear is difficult to 
approach directly. If possible, the system is simplified by 
a reduction of dimension or by a linearization technique 
(e.g., Wiggins, 1990). The linearization approach is pursued 
in the present study. The determinant of coefficients of the 
linearized system is examined to determine the first critical 
Rayleigh number (e.g., Gelaro, 1987) for both the conducting 
and insulated side wall cases. 


5.1 Local Bifurcation Theory 


Local bifurcations are those changes in stability that 
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occur for relatively simple flow patterns and solutions (e,g . , 
Argyris et al. , 1991) . For Rayleigh-B4nard convection, the 
initial conductive state and the laminar and periodic 
convective flow patterns are either steady or repetitive in 
nature. The solutions of the systems of equations derived in 
the last chapter reflect this simplicity. Procedures for 
determining the stability of these solutions are fairly 
straightforward (e.g., Wiggins, 1990). On the other hand, 
flow regimes past the periodic are more complex spatially and 
temporally. Determination of these global bifurcations is 
therefore more difficult (e.g., Wiggins, 1990). 

In the present study, for example, consider the solution 
vector (V' 54 (t) ,T 54 (t) ,ro 3 (t) ) for the case of the insulated 
side wall* For the conductive solution, all convective 
perturbations are zero, so that the solution is the zero 
vector (^ 54 (t) ,T 54 (t) ,2’o3('t) ) - (0,0,0) . This is a single 
point in the solution plot, it is time-invariant, and it is 
decidedly "local" in nature. Similarly, for fixed Rayleigh 
number, the plot of the time-invariant solution for laminar 
convection consists of two distinct points (e.g., Thompson and 
Stewart., 1986). A solution need not be time- invariant to be 
local. The periodic solution is also local because its three- 
dimensional solution plot or "trajectory" is fixed in the 
solution space and does not vary beyond one period. The 
solutions are "local" in the sense that their trajectories 
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have a decidedly simple structure in the solution plot. The 
plots of '•global'* solutions, on the other hand, have a more 
complex structure and occupy a larger portion of the solution 
space (e.g., Argyris, 1991). The loss of stability— the 
bifurcation—is the transition from one solution or structure 
to another. It is the stability of the conductive solution 
and structure in the vicinity of the first critical Rayleigh 
number that is of interest in the present study. To 
investigate the stability of the conductive regime of the 
liquid melt, the stable conductive solution must first be 
found. It is the time-invariance of this "stationary" 
solution that allows it to be determined easily. 


5.2 Stationary Solutions 

Stationary solutions are the time-invariant solutions of 
a system of ordinary differential equations. In the present 
study they are necessary for the local stability analysis 
because they are the base state in the perturbation procedure 
of the linear approximation. The definition and general 
procedure for their determination are as follows (e.g., 
Wiggins, 1990). Consider a system of ordinary differential 
equations 
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x 6H“, n€{l, 2,3, ... } (5-1) 

"Stationary" solutions are those solutions 

x€R« (5.2) 

such that 

f{x)=0 (5.3) 

stationary solutions are solutions that do not change over 
time. Each system is characterized by a set of time-invariant 
or "stationary" convective solutions representing the 
conductive and laminar convective flow regimes. This 
definition does not imply that the flow itself is physically 
stationary, but that the flow pattern is stable and does not 
vary over its own temporal scale. For example, laminar flow 
is stationary but is certainly not quiescent. The stationary 
flows represented by such stationary solutions often exhibit 
variance in the flow near flow transitions if sufficiently 
perturbed, but these flows are transient (e.g., Thompson and 
Stewart., 1986). 

In the present study, the system is expressed in terms of 
the temporal convective perturbations of the flow, , and 

of the temperature, r 3 ^ (t) or ^^^(t), and TQ^it) . The 
stationary solutions of the nonlinear system for the 
conducting side wall. 
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^«(t) = -lZ||ipr«^5^(t) + ^PrRarj^Ct) 


(5.4) 


TU(t) = -l||^r3,(t) - ^-fsAt)T^^{t) 


(5.5) 


ro3(t) = -ro3(t) + 


777 

2048 


^5<({t)ro3(t) 


(5.6) 


may be found by letting the time derivatives be zero. 


^54 <t) = r34( t) =ro3 (t) “0 (5.7) 

letting the stream function and temperature variables be 
constant , 


1^51 (t) = tit, T,tit) = Til, and T^^it) - li. 


:o3 C5.3) 

and solving for these stationary variables. The three 
solutions that result are the conductive solution, 


^4 = rf4 = 0 , To^3 = 0 
and two convective solutions. 


(5.9) 


^4 = 


rpS 

•^34 


rpS 

-^03 


1/-1692188606464 + 587059200Ra 


3^2456 596 87 




I 1489 \ 

\ 8084167 I 


(-1692188606464 + 587059200Ra) 


15 Ra 

-1692188606464 + 587059200Ra 
3287531520Ra 


(5.10) 

(5.11) 

(5.12) 


and 



(5.13) 



-1/-1692188606464 + 587 059200Ra 
3V245659687 




1489 \ 

8084167/ 


(-1692188606464 + 587059200Ra) 
isla 


_ -1692188606464 + 587059200Ra 
” ” 3 2 87 53 15 2 ORa 


(5.14) 

(5.15) 


Note that these solutions are solutions of the nonlinear 
system. The stationary conductive solution is essentially a 
zero-convection stationary solution. The two convective 
solutions differ only in the change in sign of the ^§4 and r |4 
solutions. These solutions are the two different rotational 
directions the flow may follow— flow rising or falling in the 
center of the cylinder. It should be noted that the radicals 
in these convective solutions are imaginary for a value of the 
Rayleigh number less than -2882.5. Imaginary values for the 
stream function and temperature variables indicate flows that 
do not exist. It will be shown later that this Rayleigh 
number is the critical Rayleigh number marking the transition 
between conductive and laminar convective flow. It may also 
be noted that the Prandtl number is not a parameter in these 
radicands, suggesting the possible Prandtl number- independence 
of this flow regime transition. 

For the insulated side wall case, the solutions are the 
conductive state. 
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= 0, = 0, Tti^ = 0 (5.16) 

and the two convective solutions. 
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^54 




■ ) ( -21555414014016 + 92449157 50 Ra) 


( 146286635/ 


52115 
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^34 " 

rpS.X ^ 
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4608-^ (-21555414014016 + 92449 l5750Ra) 


483875Ra 

8198031 (-21555414014016 + 9244915750Ra) 
421592157461282500Ra 
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(5.19) 


and 
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146286635/ 


4 83 87 5 Ra 


8198031 (-21555414014016 + 9244915750Ra) 
421592157461282500Ra 


(5.21) 

(5.22) 


In this case, the radicals in these convective solutions are 
imaginary for a value of the Rayleigh number less than 
-2331.6. 

In the following linear approximation of the local 
stability analysis, the linearized system of equations is 
perturbed about these stationary solutions and the first 



89 


critical Rayleigh number is determined. 


5.3 The Linear Approximation 

The two systems of ordinary differential equations that 
were derived in the last chapter from the perturbed Boussinesq 
system are strongly nonlinear. The local stability analysis 
of nonlinear systems of ordinary differential equations such 
as these is difficult to approach directly. in the present 
study, the system is simplified by a linearization technique 
(e.g., Wiggins, 1990). 

The linear approximation of the nonlinear local stability 
of the conductive solution is approached by a perturbation 
method. In this second perturbation of the original 
Boussinesq system, the basic state is the stationary 
conductive solution. The system is transformed into a system 
in terms of the convective perturbations. By using this time- 
invariant basic State, the temporal decay or growth of the 
perturbations alone indicates respectively the persistence of 
conduction or the initiation of laminar convection. 

In the linearized system of ordinary differential 
equations representing the flow, solutions are exponential 
expressions of the form 
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Ce^S C€R, X€C (5.23) 

in which X is an eigenvalue of the matrix of coefficients of 
the linearized system (e.g. , Gelaro, 1987). Each eigenvalue 
X represents one solution of the system. The sign of the real 
part of the eigenvalue determines the existence of the growth 
or the decay of its particular solution, and the magnitude of 
the real part determines the rate of this growth or decay. An 
eigenvalue with negative real part indicates exponential 
decay, a positive real part results in exponential growth, and 
a zero real part indicates neutral stability, neither growth 
nor decay. It is sufficient for any one real part of the set 
of eigenvalues to be positive for the system itself to be 
unstable (e.g., Gelaro, 1987). In the linear approximation, 
physical stability corresponds mathematically to the entire 
set of eigenvalues having negative real parts. Physical 
instability is initiated mathematically when any one of the 
real parts of the eigenvalues becomes positive. Zero real 
values of the eigenvalue therefore may reveal the changes of 
sign that in turn indicate changes in stability. 

For values of time t near zero, solutions are said to be 
"near" the closest time- invariant stationary solution. For 
the linear solutions, the perturbations ce^^ themselves are 
close to zero, and the solution approaches the stationary 
basic state. In this same neighborhood, the nonlinearities in 
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a nonlinear solution are products of very small quantities, 
effectively making them insignificant. In the vicinity of a 
stationary solution, the true nonlinear solution and the 
linear approximation are therefore asymptotic (e.g. , Wiggins, 
1990) . The result is that the linear behavior of a solution 
is the same as the nonlinear behavior near a stationary 
solution. If a linear solution is stable or unstable near the 
stationary solution, then the nonlinear solution is 
correspondingly stable or unstable (e.g., Wiggins, 1990). 

The real parts of the eigenvalues of the linear system 
consequently indicate the nonlinear local stability of the 
system. The presence of a single positive eigenvalue 
indicates linear and nonlinear instability. However, a 
problem arises for zero eigenvalues. A zero eigenvalue, which 
indicates neutral stability of the linear system, is 
indeterminate for the stability of the nonlinear system (e.g. , 
Wiggins, 1990) . Slight differences between the linear and 
nonlinear solutions near a stationary solution could "push” 
the stability of this particular solution either way. Only if 
the real parts of all eigenvalues of the linearized system are 
nonzero is the behavior of the linear approximation near a 
stationary solution the same as the nonlinear system (e.g., 
Wiggins, 1990) . 

In the physical system of the liquid melt, as well as in 
the nonlinear system modeling the flow and in its linear 
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approximation, the Rayleigh niimber is the critical parameter 
determining the flow regime stability (e.g., Kr ishnamurti , 
1973) . The Rayleigh number appears as a parameter in the 
eigenvalue equation, and it determines the sign of some of the 
eigenvalues. Critical Rayleigh numbers indicate the loss of 
stability of a stable solution, and they occur in sets of 
eigenvalues whose real parts are all negative, as the sign of 
any one real part changes to positive (e.g. , Wiggins, 1990) . 
This change in sign indicates a change in the stability of the 
system and a transition in the flow regimes. Although the 
real part of the eigenvalue is zero at this point and the 
linear approximation is not valid there, the critical Rayleigh 
number is the limit of stable and unstable solutions on either 
side that are valid. As the limit of such behaviors, it does 
indicate a change in stability. 

In the linear perturbation of the system of ordinary 
differential equations for the conducting side wall, the basic 
state is convective flow, and the perturbation represents a 
small disturbance of this flow. In general, the basic state 
is a stationary flow or solution, and the perturbation is a 
function of time (e.g., Gelaro, 1987): 

^S4(t) = ^54 *■ (5.24) 


T34(t ) = 2/4 + 2^4 (t) 


( 5 . 25 ) 
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To3 ( t ) * r /5 + Tijit) (5.26) 

These are substituted into the nonlinear system, 

^54(t) = - J:Z| | Q pr^ 5 ^(t) + i|prRar 3 ^(t) (5.27) 



ifiCt) = fLit) - - M|Z|ir/4(t) - Mtlrl^Ti^it) 

(5.31) 

T«(t) = ^^r/ 4 ^l 4 (t) * - 48ro^(t) (5.32) 

Note that the system is now expressed in the perturbation 
variables. The stability of these variables determines the 
stability of the entire system. It is the eigenvalues of 
this system that are to be investigated. It is more 
illustrative of the method to determine these by transforming 
the system into a system of variational equations than to find 
the eigenvalues directly (e.g., Gelaro, 1987). The system 
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above has perturbation solutions of the form; 


^Lit) 


(5.34) 

Tiiit) 


(5.34) 

Tq3 (t ) 

= ^ 
*‘ 03 ^ 

(5.35) 


The effect of the sign of the eigenvalue on the growth of the 
perturbation is evident here. If these expressions are 
substituted into the linearized system, the variational 
equations become; 

+ •||prRaf34 = 0 (5.36) 
(l - - (~|||^ + A)f34 - = 0 (5.37) 

- <■‘8 + A)f„3 = 0 (5.38) 

The system above allows determination of the eigenvalues for 
any stationary solution. In the present study the conductive 
stationary solution is of interest, and the system can be 
simplified considerably by letting the stationary solutions be 
the zero-convection solution; 
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-If = 0 

(5.39) 

^ _ ( 138728 ^ _ 

\ 4095 ° 

(5.40) 

(48 + X) 2^03 = 0 

(5.41) 


For this system to have a nontrivial solution, the determinant 
of coefficients must be zero. The eigenvalue equation is a 
third degree polynomial in X: 


(48 + X)[A ,2 + i/lZ 868 Pr ^ 138728 \ ^ 17868Pr 138728 
[ V 355 4095 j 355 4095 


42 

71 


PrRa 


= 0 


The eigenvalue solutions of this equation are: 


= -48 


(5.42) 


(5.43) 


1/ 17 868 Pr ^ 138728 \ If/ 17 868 Pr ^ 138728 \2 
2 \ 355 4095 j 2 [\ 355 4095 / 


7 17868 138728 
\ 355 4095 




(5.44) 


Eigenvalue X^^ is a constant, negative real number, indicating 
stability for its particular solution. Because the Prandtl 
mamber is always positive, the real part X 3 will always be 
negative, the radical being either a positive real number or 
a pure imaginary. 
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The sign of the real part of X 2 , however, depends upon 
the radicand. Because X 2 is of the form 

-fi + - fg (5.45) 


the sign of the quantity 


7 17 868 138728 
[ 355 4095 



(5.46) 


determines the sign of the real part of X 2 and consequently 
the stability of the system. if this quantity is positive, 
which occurs for a Rayleigh number less than ~2882.5, the 
radical term will be either a positive real number less than 


^ 1/ 17868PI ^ 1387 28 \ ,5 
2 \ 355 4095 / ( ' ) 

or an imaginary number, and the real part of X 2 will be 
negative. In this case, all the eigenvalues are negative and 
the system is stable. If the quantity is negative, for a 
Rayleigh number greater than ~2882.5, then the radicand will 
be a positive number such that the real part of X 2 will be 
positive. The system is then unstable. This value of -2882.5 
is the first critical Rayleigh number for the conducting side 
wall case. The condition that the critical Rayleigh number be 
at the change in sign of the real part of the eigenvalue has 
been met. The supposition in section 5.2 concerning the value 
of the critical Rayleigh number based on the stationary 
solutions was correct. 



97 


It is evident that the Prandtl number has no effect on 
the sign of X 2 or changes in the stability of the system. 
Therefore the change in stability is a function of the shape 
of the container and not of the fluid properties. In this 
linear model of stability, the Prandtl number does however 
have an effect on the rate of growth of the convective 
perturbation near the instability. The higher the Prandtl 
number, the larger the positive real eigenvalue and the 
greater the growth of the perturbation. 

The behavior of the system with respect to the transition 
between stationary solutions may be seen easily in the 
bifurcation diagrams of the stationary solutions for the 
stream function and temperature temporal coefficients (Figs. 
22, 23, and 24). 

The procedure for determining the critical Rayleigh 
number is repeated for the case of insulated sides (see 
Appendix A for the calculations) . After considerable algebra, 
the critical Rayleigh number in this case is found to be 
~2331.6. The bifurcation diagrams are much the same as those 
for the conducting side wall case (Figs. 25, 26, and 27). 
This lower value of the critical Rayleigh number validates the 
model's treatment of the thermal boundary conditions. From a 
physical perspective, the lower value for the instability here 
indicates that less of a temperature gradient is required to 
destabilize the fluid colximn. This is exactly the case for 
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the Insulated side wall, through which there is no heat loss. 
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Figure 23. Bifurcation diagram for r 34 (t) , conducting side 


wall . 
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Ficnire 24. Bifurcation diagram for Tq 3 (t) , conducting side 


wall. 
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gicmre 25 . Bifurcation diagram for 1^54 (t), insulated side 


wall. 
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ire 26 . Bifurcation diagram for 1*54 (t), insulated side 



2331.6 

RAYLEIGH NUMBER, Ra 


Bifurcation diagram for Tn^it) , insulated side 






CHAPTER 6. DISGPSSION OF RESULTS 


The Chebyshev-GalerJcin spectral model of Rayleigh-B4nard 
convection in this study is formulated upon, and exhibits, 
certain features that have been observed in experiments. 
Axisymmetric flow is found in cylinders of equal height and 
radius (e.g., Muller et al., 1984). The convergence of this 
type of flow at the center of the cylinder moves the center of 
rotation of the single convection cell outward past R/2, the 
midpoint of the radius (e.g., Mtiller et al., 1984). Although 
the behavior of a low-Prandtl nximber fluid is investigated in 
this study, both experiment and the model itself demonstrate 
the Prandtl number independence of the first critical Rayleigh 
number (Krishnamurti , 1973) . 

The results for the first critical Rayleigh mamber from 
this study are both slightly higher than those found by 
Char Ison and Sani (1970) in their Rayleigh-Ritz formulation 
for the same problem. For the conducting side wall, the value 
of 2882.5 is 13.3% higher than their value of 2545.0, and for 
the insulated side wall, the value of 2331.6 is 3.1% higher 
than their value of 2261.9. 

The experimental value of the first critical Rayleigh 
number for liquid gallium determined by Miiller et al. (1984) 
had to be interpolated from their graph (Fig. 28) for their 
aspect ratio of h/d=0.5 (y=l) because the authors did not 
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state the exact numerical result. Interpolation from the 
graph gives a value for of approximately 2700 for the 
Bridgman configuration (insulated side wall) . Although the 
value of 2331.6 from the present study differs by 13.6% from 
this experiment, it is closer than that found by Charlson and 
Sani (1970) . 



Rayleighnumber Ra 
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Ficnire 28. Graph from Mtiller et al. (1984) showing the flow 
regimes and transitions for liquid gallium. 



CHAPTER 7. CONCLUSION 


In this study, the onset of laminar axisymmetric 
Rayleigh-B^nard convection is investigated for a fluid in a 
cylindrical container whose radius is equal to its height. A 
simplified model of the vertical Bridgman crystal growth 
configuration is used. Two different cases are considered for 
the thermal boundary conditions at the side wall; conducting 
and insulated surfaces. All surfaces are assumed to be solid 
and no~slip. 

The governing Boussinesg system is first perturbed. It 
is then simplified by introducing a Stokes stream function and 
by taking the curl of the Navier*Stokes equations. A 
Chebyshev-Galerkin spectral model reduces the simplified 
system to a system of first-order ordinary differential 
equations. A local stability analysis using the linearized 
system determines the two values of the first critical 
Rayleigh number for the insulated and conducting side walls. 

Although this study investigates the onset of laminar 
convection for a low-Prandtl number liquid metal in a cylinder 
with equal radius and height, the results are more general. 
The critical Rayleigh number at this transition is independent 
of the Prandlt niamber, so that the Ragi values obtained here 
are valid for all fluids within the parameters of the 
formulation. Furthermore, the presence of a single convection 


108 



109 


cell in cylinders of aspect ratio higher than the value of one 
in this study (e.g. , Tritton, 1988) suggests that the model is 
applicable to the determination of RSci for these high-aspect 
ratio cylinders. 

The relative magnitudes of the first critical Rayleigh 
ntimber values that are obtained, 2882.5 for the conducting 
side wall and 2331.6 for the insulated case, tend to validate 
the method. The Ra^i for the conducting side wall should be 
higher than that for the insulated side wall, the heat lost 
through the conducting surface not being available to initiate 
the flow transition. 

The approach used in this study results in a value of 
Ra^^ for the case of the insulated side wall that is closer 
than that of the niimerical work of Charlson and Sani (1970) to 
the experimental value determined by Muller et al. (1984) . 

A more precise approach to the problem presented in this 
study would be to determine the flow and temperature 
distributions in the laminar convective regime by some 
numerical method rather than by curve-fitting. These 
distributions could be approximated by polynomials and 
Chebyshev series . The rest of the procedure would be the 
same. 

The procedure used in Chapter 5 to determine the first 
critical Rayleigh number may also be used to determine the 
second critical Rayleigh number. In this case the nonzero 
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laminar convective stationarY solutions are used in the 
linearized system and its eigenvalue equation. Although the 
eigenvalue equation is a numerically more complicated cubic, 
its solution is still possible (e.g., Gelaro, 1987). 

The Chebyshev-Galerkin spectral method used in this study 
to determine the first critical Rayleigh number for a low- 
Prandtl number crystal melt must be modified to find the full 
sequence of critical Rayleigh numbers. The Chebyshev series 
expressions for the stream function and the temperature need 
to be considerably longer to account for the increased number 
of flow modes or length scales that need to be represented as 
the flow become more turbulent. These expressions must also 
allow inclusion of the inertial term of the Navier-Stokes 
equation to account for the strong nonlinearity for any non- 
zero velocity field. A serious problem with the extrapolation 
of the specific method used here is the difficulty of 
satisfying the boundary conditions. The coefficients of the 
spatial distributions in the present study could be chosen to 
satisfy the boundary conditions because there was only one 
simple flow field to model. If these coefficients are 
unknown, the boundary conditions are not automatically 
satisfied and impose added constraints. 



REFERENCES 


Argyris, J., Faust, G., and Haase, M. (1991). "Xaog, an 
adventtire in chaos . '• Computer Methods in Applied Mechanics and 
Engineering 91, 997-1091. 

B^nard, H. (1900) . "Les tourbillons cellulalres dans une nappe 
liquide." Revue G^ndrale des Sciences Pares et Appliqu6es ll, 
1261-1271. 

B§nard, H. (1901) . ”Les tourbillons cellulaires dans une nappe 
liquide transportant de la chaleur par convection en regime 
permanent.” Annales de Chimie et de Physique 23, 62-144. 

Berge, P. , Pomeau, Y. , and Vidal, C. (1984) . Order within 
Chaos. John Wiley & Sons, Inc., New York, N.Y. 

Block, M. J. (1956) . "Surface tension as the cause of Benard 
cells and surface deformation in a liquid film.” Nature 178, 
650-651. 

Boussinesq, J. (1903) . Th^orie Analytique de la Chaleur, Vol. 
2. Gauthier-Villars, Paris. 

Canute, C. , Hussaini, M.Y., Quarteroni, A., and Zang, T.A. 
(1990). Spectral Methods in Fluid Mechanics. Springer-Verlag, 
New York, N.Y. 

Carlson, F.M., Fripp, A.L., and Crouch, R.K. (1984). "Thermal 
convection during Bridgman crystal growth. ” Journal of Crystal 
Growth 68, 747-756. 

Chandrasekhar, S. (1961). Hydrodynamic and hydromagnetic 
stability . Clarendon Press, Oxford. 

Charlson, G.S. and Sani, R.L. (1970). "Thermoconvective 
instability in a bounded cylindrical fluid layer." 
International Journal of Heat and Mass Transfer 13 , 1479-1496. 

Crouch, R.K., Fripp, A.L., Debnam, W.J. , Clark, I.O., Barber, 
P.G., and Carlson, F.M. (1985). "Experimental investigation of 
the effects of gravity on thermosolutal convection and 
compositional homogeneity in Bridgman grown, compound 
semiconductors." Acta J^tronautica 12, 923-929. 

Davis, S.H. (1968). "Convection in a box: on the dependence of 
preferred wave-number upon the Rayleigh number at finite 
amplitude.” Journal of Fluid Mechanics 32, 619-624. 


Ill 



112 


Deane, A.E. and Sirovich, L. (1991) . "A computational study of 
Rayleigh-B4nard convection. Part 1. Rayleigh-number Scaling." 
Journal of Fluid Mechanics 222, 231-250. 

Drazin, P.G. and Reid, W.H. (1981) . Hydrodynamic Stability. 
Cambridge University Press, Cambridge. 

Gelaro, R. (1987). "Linear stability analysis," in: Nonlinear 
Hydrodynamic Modeling: A Mathematical Introduction , Ed. H.N. 
Shirer, Springer-Verlag, New York, N.Y. 

Gleick, J. (1988) . Chaos: Making a New Science. Penguin Books, 
New York, N.Y. 

Haldenwang, P. (1986) . "Unsteady numerical simulation by 
Chebyshev spectral methods of natural convection at high 
Rayleigh number." Proceedings of the 1986 ASME Winter Annual 
Meeting, "Significant Questions in Buoyancy-Affected 
Enclosures or Cavity Flows," Anaheim, CA, 45-51. 

Higgins, R.W. (1987). "From the equations of motion to 
spectral models," in: Nonlinear Hydrodynamic Modeling: A 
Mathematical Introduction, Ed. H.N. Shirer, Springer-Verlag, 
New York, N.Y. 

Horsch, G.M. (1988) . "Cooling- induced convective littoral 
currents in lakes," dissertation presented to the University 
of Minnesota, in partial fulfillment of the requirements for 
the Degree of Doctor of Philosophy. 

Kim, K.M. , Witt, A.F. , and Gatos, H.C. (1972). "Crystal growth 
from the melt under destabilizing thermal gradients." J. 
Electrochem. Soc, 119, 1218. 

Knuteson, D. J. (1989) . "Experimental study of convection in 
tin in a vertical Bridgman configuration," dissertation 
presented to the Department of Chemical Engineering, 
University of Florida, in partial fulfillment of the 
requirements for the Degree of Doctor of Philosophy. 

Koschmieder, E. (1966) . "On convection on a uniformly heated 
rotating plane." Beitr. Phys. Atmos. 39, 1-11. 

Koschmieder, E. (1967). "On convection on a uniformly heated 
rotating plane." Beitr. Phys. Atmos. 40, 216-225. 

Koschmieder, E. (1969) . "On the wavelength of convective 
motions." J. Fluid Mech. 350, 527-530. 



113 


Kr ishnamurti , R. (1970a). "On the transition to turbulent 
convection. Part 1.” Journal of Fluid Mechanics 42, 295-307. 

Kr ishnamurti, R. (197 Ob). "On the transition to turbulent 
convection. Part 2.” Journal of Fluid Mechanics 42, 309-320. 

Krishnamurti, R. (1973), "Some further studies on the 
transition to turbulent convection." Journal of Fluid 
Mechanics 60, 285-303. 

Le Queri, P. and Alziary de Roquefort, T. (1985) . "Computation 
of natural Convection in two-dimensional cavities with 
Chebyshev polynomials," Journal of Computational Physics 57, 
210-228. 

Liang, S.F., Vidal, A., and Acrivos, A, (1969). "Buoyancy- 
driven convection in cylindrical geometries." Journal of Fluid 
Mechanics 32, 619-624. 

Lorenz, E.N. (1963), "Deterministic nonperiodic flow." J, 
Atmos. Sci. 20, 130-141. 

Miiller, G. , Neumann, G. , and Weber, W. (1984). "Natural 
convection in vertical Bridgman configurations". Journal of 
Crystal Growth 70, 78-93. 

Mullins, W.W. and Sekarka, R.F, (1964). "Stability of a planar 
interface during solidification of a dilute binary alloy." J. 
Appl. Phys. 35, 444. 

Nese, J.N. (1987). "The transition to turbulence," in: 
Nonlinear Hydrodynamic Modeling: A Mathematical Introduction, 
Ed. H.N. Shirer, Springer-Verlag, New York, N.Y. 

Orszag, S.A. (1971a) . "Numerical simulation of incompressible 
flows within simple boundaries. I. Galerkin (spectral) 
representations." Studies in Applied Mathematics L(4) , 293- 
327. 

Orszag, S.A. (1971b) . "Galerkin approximations to flows within 
slabs, spheres, and cylinders." Physical Review Letters 
26(18), 1100-1103. 

Ostrach, S. (1985). "Fluid mechanics in crystal growth— The 
1982 Freeman Scholar Lecture", Journal of Fluids Engineering 
105, 5-20. 

Ostrach, s., and Pneuli, D. (1963). "The thermal instability 
of completely confined fluids inside some particular 



114 


configurations” Journal of Heat Transfer 85, 1346-1354. 

Parker, S.G. and Johnson, R.E. (1981) in: Preparation and 
Properties of Solid State Materials , Eds. W.R. Wilcox and R.A. 
Levever, Marcel Oekker, New York, N.Y. 

Pearson, J.R.A. (1958) . On convective cells induced by surface 
tension." Journal of Fluid Mechanics 4, 489-500. 

Pellew, A. and Southwell, R.V. (1940). "On maintained 
convective motion in a fluid heated from below." Proa. R. Soc. 
176A, 312-343. 

Lord Rayleigh (1916) . "On convective currents in a horizontal 
layer of fluid when the higher temperature is on the under 
side." Phil. Mag. 32, 529-546, 

Saltzmann, B. (1962) . "Finite amplitude free convection as an 
initial value problem— I." J. Atmos. Sci. 19, 329-341. 

Segel, L.A. (1969) . "Distant side-walls cause slow-amplitude 
modulation of cellular conection." Journal of Fluid Mechanics 
38, 293^324. 

Shirer, H.N. (1987). "A simple nonlinear model of convection," 
in: Nonlinear Hydrodynamic Modeling: A Mathematical 
Introduction, Ed. H.N. Shirer, Springer-Verlag, New York, N.Y. 

Tarman, I.H. (1989). "Analysis of turbulent thermal 
convection," dissertation presented to the Department of 
Applied Mathematics, Brown University, in partial fulfillment 
of the requirements for the Degree of Doctor of Philosophy. 

Thompson, J.M.T. and Stewart, H.B. (1986) , Nonlinear Dynamics 
and Chaos. John Wiley & Sons Ltd., Chichester. 

Tiller, W.A., Jackson, K.A., Rutter, J.W. , and Chalmers, B. 
(1953) . "The redistribution of solute atoms during the 
solidification of metals." Acta Metallurgica 1, 428. 

Tritton, D. J. (1988) . Physical Fluid Dynamics, 2nd ed. 
Clarendon Press , Oxford . 

Wiggins, S. (1990) . Introduction to Applied Nonlinear Systems 
and Chaos. Springer-Verlag, New York, N.Y. 

Zierep, J. (1963). "Zur Theorie der Zellularkonvektion V." 
Beitr. Phys. Atmos. 36, 70-76. 



APPENDIX A. LOCAL STABILITY ANALYSIS. INSUIATED SIDE WALL 

The local stability analysis of the system of nonlinear 
ordinary differential equations obtained in Chapter 4 for the 
case of the insulated side wall is similar to that for the 
conducting side wall. The system is perturbed about the 
conductive stationary solution and is linearized. The first 
critical Rayleigh number is determined from the eigenvalue 
equation of the linearized system. The procedure follows that 
of Gelaro (1987). 

In the determination of the local stability of a 
stationary convective solution, the basic state is the 
stationary convective flow, and the temporal perturbation 
represents a small disturbance of this flow. As in the case 
of the conducting side wall, the independent variables for the 
insulated side wall are represented by first-order expansions: 


= ^54 + ^54 it.) 

(A.l) 

Tsiit) = r /4 + Tliit) 

(A.2) 

Tos(t) = Tis + Tfsit) 

(A. 3) 


These expressions are substituted into the nonlinear system 
derived for the case of the insulated side wall (eqns. 4.92- 
4.94), 
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Vs4i^) = - 


17 86 8 pj. ^ j. j ^ pr RaT54 ( t ) 


355 


(A. 4) 


r’4(t) = 


114636 ^ 38915152 „ 

TSeHor^^^ 


1205993 


3209800 

6029965 


^54(t)ro3(t) 


(A.5) 


r^3(t) =-48To3(t) * |||2|^,^(t)r,^(t) 


(A. 6) 


The resulting system is linearized by retaining only those 
terms of the order of the perturbations: 


i^4(t) - - 


17868 

355 


Pr^4(t) + 3^PrRaT/4<t) 


426 


(A. 7) 


T^P' - 114636 /t-N - 3209800 rpS 

■ 1205993’'"'"’ 6029963^'’=’^"'"' 


38915152 mP 
6029965 


3209800 .,S rnP 


6029965 




(A. 8) 




16384 


^ lllff -48To^(t) 


(A. 9) 


This system of first-order ordinary differential equations has 
exponential solutions of the form: 




(A. 10) 


Ti^ (t) 


= T 

•^54® 


(A. 11) 


r/aCt) = foae^t 


(A. 12) 
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Substitution of these expressions into the linearized system 
yields the variational system of equations: 

.(l|||8pj. , ^ - 0 (A. 13) 


114636 _ 3209800 p _ / 38915152 ^ 

1205993 6029965 ^ 1361605 ) ‘ 


3209800 .s A _o 
6029965 ° 


(A.14) 


63805 

16384 



+ 


63805 

16384 


- 


(48 + A.) fo3 * 0 


(A.15) 


The system above allows determination of the eigenvalues for 
any stationary solution. In the present study, the stability 
of the conductive stationary solution for the insulated side 
wall. 


^4 = 0, ifi = 0, ifa = 0 (A. 16) 

is of interest, and the system can be simplified by letting 
the stationary solutions be this zero-convection solution: 




(A. 17) 


114636 ^ _ p8915152 

1205993 ^®^ \ 1361605 



(A. 18) 


(48 + X) fo3 =* 0 (A. 19) 

For this system to have a nontrivial solution, the determinant 
of coefficients must be zero. This eigenvalue equation is a 
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third degree polynomial in X: 


(48 + X) 


X^ + 


17868PI 

355 


^ 38915152 \ ^ 17868Pr 38915152 
1361605 / 355 1361605 


/ 2765 y 114636 
V 426 A 1205993 


1 

PrRa 

/ 


0 


(A. 20) 


The signs of the real parts of the eigenvalue solutions 
of this equation determine the stability of the system. If 
the real parts of all the eigenvalues are negative, then all 
linear perturbations decay, and the system is stable. If the 
real part of any one eigenvalue is positive, its particular 
perturbation solution grows without bound, making the whole 
system unstable. The transition from stability to instability 
occurs at the change in sign of a single real part in a set of 
negative real parts. The first critical Rayleigh number, 
which marks this transition from stability to instability, is 
found by varying the Ra parameter to obtain the first zero 
real part at which the sign changes. 

The eigenvalue solutions of the above equation are; 


= -48 


(A. 21) 


1/ 17868Pr ^ 38915152 \ If/ 17 868Pr ^ 38915152 

2\ 355 1361605 ] 2[\ 355 1361605 ] 


4Pr 


/ 17868 

38915152 

2765 

114636 p \1 

\ 355 

1361605 

426 

1205993 1 


(A. 22) 


Eigenvalue is a constant negative real number, indicating 
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stability for its particular solution. Because the Prandtl 
number is always positive, the real part of X 3 will always be 
negative, the radical being either a positive real niamber or 
a pure imaginary. 

The sign of the real part of X 2 , however, depends upon 
the radicand. Because X 2 is of the form 

-fi + - fg (A. 23) 


the sign of the quantity 


4Pr 


17 86 8 y 38915152 \ 
355 j\ 1361605 / 


2765 y 114636 W ' 
426 A 1205993/ 


(A. 24) 


determines the sign of the real part of X 2 and consequently 
the stability of the system. If this quantity is positive, 
which occurs for a Rayleigh number less than ~2331.6, the 
radical term will be either a positive real number less than 


II 17868Pr ^ 3 8915152 \ 
2\ 355 1361605 / 


(A.25) 


or an imaginary number, and the real part of X 2 will be 
negative. In this case, all the eigenvalues are negative and 
the system is stable. If the quantity is negative, for a 
Rayleigh number greater than -2331.6, then the radicand will 
be a positive number such that the real part of X 2 will be 
positive. The system is then unstable. This value of -2331.6 
is the first critical Rayleigh momber for the insulated side 
wall case. The condition that the critical Rayleigh number be 
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at the change in sign of the real part of the eigenvalue has 
been met. The supposition in section 5.2 concerning the value 
of the critical Rayleigh number based on the stationary 
solutions (egs. 5.17-’5.19 and 5.20-5.22} has again been shown 
to be correct. In this case also, the first critical Rayleigh 
nximber is independent of the Prandtl number. 



